#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "IDoubleMatrix.h"

namespace LatoolNet {
	using namespace System;

	ref class DoubleSymmetricTridiagonalMatrix : public IDoubleMatrix {
	private:
		double * d;
		double * e;
		long * ipiv;
		int m_rownum;
		int m_colnum;
		bool m_isFactorized;
		bool m_isPositiveDefinite;

		DoubleSymmetricTridiagonalMatrix(DoubleSymmetricTridiagonalMatrix ^ orig) {
			m_rownum = orig->m_rownum;
			m_colnum = orig->m_colnum;

			d = new double[orig->m_rownum];
			e = new double[orig->m_rownum - 1];

			for (int i = 0; i < orig->m_rownum; i++) {
				d[i] = orig->d[i];
			}
			for (int i = 0; i < orig->m_rownum - 1; i++) {
				e[i] = orig->e[i];
			}
			ipiv = new long[m_rownum];

			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}

			m_isFactorized = orig->m_isFactorized;

		};
	public:
		DoubleSymmetricTridiagonalMatrix(int rownum, int colnum) {
			if (rownum != colnum) {
				throw gcnew LatoolException("Number of rows and columns must be equal for symmetric matrix.");
			}
			m_rownum = rownum;
			m_colnum = colnum;

			d = new double[rownum];
			e = new double[rownum -1];

			for (int i = 0; i < rownum; i++) {
				d[i] = 0.0;
			}
			for (int i = 0; i < rownum - 1; i++) {
				e[i] = 0.0;
			}
      ipiv = new long[rownum];

			m_isFactorized = false;
		};

		~DoubleSymmetricTridiagonalMatrix(){
			delete[] d;
			delete[] e;
		};

		virtual IDoubleMatrix ^ Clone() {		
			return (IDoubleMatrix^) gcnew DoubleSymmetricTridiagonalMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property double default[int, int] {
			double get(int i, int j) {
				if (i == j) {
					return d[i];
				} else if (i == j + 1) {
					return e[j];
				} else if (j == i + 1) {
					return e[i];
				} else {
					return 0.0;
				}
			}
			void set(int i, int j, double value) {
				if (i == j) {
					d[i] = value;
				} else if (i == j + 1) {
					e[j] = value;					
				} else if (j == i + 1) {
					e[i] = value;
				} else {
					return;
				}
				m_isFactorized = false;
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::DoubleSymmetricTridiagonal; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
		}

		virtual void Invert() {
			throw gcnew LatoolException("Inversion of tridiagonal matrix is not supported by LAPACK.");
		};

		virtual void Factorize() {

			long n = m_rownum;
			long info;

			dpttrf_(&n, d, e, &info);

			if (info == 0) {
				m_isFactorized = true;
				return;
			}

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			} else {
				throw gcnew LatoolException("This matrix is not positive definite. Symmetric tridiagonal indefinite matrix is not supported by LAPACK. Please use Matrix::DoubleSymmetrix."); 
			}

		};

		virtual void Solve(IMatrix ^ b) {

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			long n = m_rownum;
			long nrhs = b->ColNum;
			long ldb = Math::Max(1, n);
			long info;

			dpttrs_(&n, &nrhs, d, e, gb->dat, &ldb, &info);					

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			long n = m_rownum;
			long nrhs = b->ColNum;
			long ldb = Math::Max(1, n);
			long info;

			dptsv_(&n, &nrhs, d, e, gb->dat, &ldb, &info);
			
			if (info == 0) {
				m_isFactorized = true;
				return;
			}

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			} else {
				throw gcnew LatoolException("This matrix is not positive definite. Symmetric tridiagonal indefinite matrix is not supported by LAPACK. Please use Matrix::DoubleSymmetrix."); 
			}
		};


	};
}